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The dynamical low-rank approximation of time-dependent matrices is a low-rank 
factorization updating technique. It leads to differential equations for factors of the 
matrices, which need to be solved numerically. We propose and analyze a fully ex- 
plicit, computationally inexpensive integrator that is based on splitting the orthogonal 
projector onto the tangent space of the low-rank manifold. As is shown by theory and 
illustrated by numerical experiments, the integrator enjoys robustness properties that 
are not shared by any standard numerical integrator. This robustness can be exploited 
to change the rank adaptively. Another application is in optimization algorithms for 
low-rank matrices where truncation back to the given low rank can be done efficiently 
by applying a step of the integrator proposed here. 

1 Introduction 

Low-rank approximation of large matrices and tensors is a basic model reduction technique in a 
wide variety of applications ranging from quantum physics to information retrieval. In the present 
paper, we consider the low-rank approximation of time-dependent matrices A(t), to < t < t, which 
are either given explicitly or are the unknown solution of a differential equation A(t) = F(A(t)). 
In the first case, instead of performing an (expensive) low-rank approximation via singular value 
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decompositions independently for every t, one would prefer a computationally inexpensive updat- 
ing procedure that works only with the increments of A and is free from decompositions of large 
matrices. In the second case, one would like to find an approximate solution to the differential 
equation working only with its low-rank approximation. 

Both cases can be dealt with by dynamical low-rank approximation where A{t) is approximated 
by Y(t) of fixed rank r by imposing that its time derivative should satisfy 

||F(f)-A(f)|| =min, (1) 

where the minimization is over all matrices that are tangent to Y(t) on the manifold .4( r of matrices 
of rank r, and the norm is the Frobenius norm. In the case of a differential equation A = F(A), the 
above minimization is replaced with 

||y(0-F(y(0)||=min. (2) 

In both cases, this leads to a differential equation on the manifold of matrices of a given rank r, 
which is then solved numerically. 

For large time-dependent matrices, this approach was proposed and analyzed in Q, and first 
numerical applications were given in 021. The approach and its error analysis were extended to 
tensors in the Tucker format in 0, Very recently, in HI El S3 the dynamical low-rank approxi- 
mation approach was extended to tensors in the tensor train (TT) format studied in [13] and the 
hierarchical Tucker (HT) format studied in 0. In quantum dynamics, such an approach was used 
previously in the multiconfiguration time-dependent Hartree (MCTDH) method IfTOl [TT1I to de- 
termine an approximate solution to the time-dependent multi-particle Schrodinger equation, and 
the above minimization process is known there as the Dirac-Frenkel time-dependent variational 
principle (see, e.g., Q). 

The dynamical low-rank approximation leads to differential equations that need to be solved 
numerically. In this paper we present a numerical integration technique that is fully explicit and, in 
contrast to standard integrators such as (explicit or implicit) Runge-Kutta methods, does not suffer 
from a possible ill-conditioning of certain matrices arising in the differential equations. This new 
method is based on a splitting of the projector onto the tangent space of the low-rank manifold 
at the current position. A different splitting algorithm was recently proposed in iTPfll . where the 
differential equations are split along different components. The projector splitting discussed here 
offers, however, several advantages: it leads to a much simpler and less expensive time- stepping 
algorithm and it can be shown to enjoy remarkable robustness properties under the ill-conditioning 
mentioned before. 

In Section[2]we recapitulate the dynamical low-rank approximation of matrices, and in Section[3] 
we describe the novel splitting integrator. In Section [4] we analyze the robustness under over- 
approximation, that is, under the dreaded ill-conditioning mentioned above. Section [5] presents 
numerical experiments. The final section addresses some perspectives for the use of the proposed 
integrator and extensions. 

2 Dynamical low-rank approximation of matrices 

Let r be a given rank, and let M r denote the manifold of all real mxn matrices of rank r (typically, 
r <C m,n). In the dynamical low-rank approximation to matrices A(t) G IR mx ", the approximation 
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Y(t) G ^ r is represented in a non-unique way as 

Y(t) = U(t)S(t)V(t) T , (3) 

where U(t) G IR mx '" and V(t) G W ixr each have r orthonormal columns, and S(t) G W xr is an 
invertible matrix. This looks similar to the singular value decomposition, but S(t) is not assumed 
diagonal. 

It is shown in [5, Prop. 2.1] that, given such a decomposition Yq = UqSqVq T of the starting value 
and imposing the gauge conditions 

U(t) T U(t)=0, V(t) T V(t) = 0, (4) 

the solution Y(t) G to the time-dependent variational principle ([I]) admits a unique decompo- 
sition ([3]), and the factors U(t),S(t),V(t) satisfy the following system of differential equations: 

U(t) = (I-U(t)U(t) T )A(t)V(t)S(t)- 1 

V(t) = (/ -V{t)V(t) T )A(t) T U{t)S{t)- J (5) 
S(t) = U(t) T A(t)V(t). 

This system of differential equations is to be solved numerically. The numerical solution of ([5]) by 
standard integrators (e.g., of Runge-Kutta type) becomes cumbersome if S(t) is nearly singular. 
This situation occurs in the case of over-approximation, when the true rank of the solution (or 
approximate rank) is smaller than the chosen rank r. This is a realistic case: the effective rank 
may not be known in advance, and it is often reasonable to overestimate the rank for accurate 
approximation. To avoid the possible singularity of S(t) usually some sort of regularization is 
used. However, such a regularization introduces errors that are poorly understood. 

The minimization condition ([!]) states that at an approximation Y{t) G the derivative Y(t) is 
obtained by orthogonally projecting A (t) onto the tangent space T Y ( t )-^r of the rank-r manifold at 
Y(t): 

Y(t)=P(Y(t))A(t), (6) 

where P(Y ) is the orthogonal projector onto the tangent space Ty^ r of the manifold at Y G ^t r . 
The projector has a simple representation [5, Lemma 4.1]: for Y = USV T as in 

P(Y)Z = ZVV T -UU T ZVV T + UU T Z. (7) 

Note that UU T is the orthogonal projector onto the range M(Y) of Y = USV T , and VV T is the 
orthogonal projector onto the range &(Y T ), so that we can also write 

P(Y)Z = ZP M( YT^—P S g^ Y ) Z P^(Y T ) +P,<%(Y) Z - (8) 

3 The integrator 

3.1 First-order splitting method, abstract formulation 

Let a rank-r approximation Yq to A(to) be given and consider a step of the standard Lie-Trotter 
splitting of <|6]) with ^ from t Q to t\ = t + h: 
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1. Solve the differential equation Yj = AP^, Y ^) w i tn initial value Yj(to) = Yq on the interval 

to < t < t\. 

• • 

2. Solve the differential equation Y u = —P^^ y „)^P^(y t ) w ^ tn initial value F//(?o) = Y/(t\) on 
the interval to < t < t\ . 

• • 

3. Solve the differential equation Y m = Pgg(y\A with initial value Yjjj(to) = Yjj(t\) on the 
interval to < t < t\ . 

Finally, take Y\ = Yni{t\) as an approximation to Y(t\), the solution of ^ at t\. By standard theory, 
this is a method of first-order accuracy. Remarkably, each of the split differential equations can be 
solved exactly in a trivial way. 

Lemma 3.1. The solution of 1 . is given by 

Y,(t)=U 1 (t)S 1 (t)V I (t) with (UjSi)' =AVi, V> = 0. (9) 

The solution of 2. is given by 

Y u (t) = U n (t)S n (t)Vn(t) with S U = -UjjAVu, U n = 0, V n = 0. (10) 

The solution of 3. is given by 

Y m (t) = U ni (t)S m (t)Vni(t) with (V ni Sj n )' = A U IIh U m = 0. (11) 

Proof. We first notice that each of the terms on the right-hand side of ([8]) is in the tangent space 
Ty^r, because for the first term the orthonormality V T V = I yields 

P(Y)(ZVV T ) =ZVV T + UU T ZVV T -UU T ZVV T =ZVV T , 

so that ZVV T E T Y Jlr- Similarly we also have UU T Z e T Y Jt r and UU T ZVV T E T Y Jt r . It 
therefore follows that the solutions of 1.-3. all stay of rank r. Hence, Yj(t) can be factorized as 
Yj(t) — Ui(t)Si(t)Vi(t) with an invertible sxs matrix 5/ and t//,V/ having orthonormal columns. 
All the matrices can be chosen to be differentiable, and we thus have 

Y 1 = (U,S 1 )'V^ + (U I S I )V^ 

• • • • 

which by 1. must equal Yj =AViV I . We observe that this is satisfied if (UiSi) =AV/and V/ = 0. 
The proofs for Steps 2. and 3. are similar. □ 



3.2 First-order splitting method, practical algorithm 

Lemma 3.1 leads us to the following algorithm. Given a factorization ([3]) of the rank-r matrix 
Yq = UqSoVq T and denoting the increment AA = A(t\) — A(to), one step of the integrator reads as 
follows: 

1. Set 

K x = U S + AAV 

and compute the factorization 

with U\ having orthonormal columns and anrxr matrix S\ (using QR or SVD). 
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2. Set 



3. Set 



and compute the factorization 



S = 5i-£/i T MV . 



Li = V S T +AA T U l 



ViSj=L h 



with V\ having orthonormal columns and anrxr matrix Si (using QR or SVD). 
The algorithm computes a factorization of the rank-r matrix 

which is taken as an approximation to Y(t\). Note that Y\ is identical to the result of the abstract 
splitting algorithm of Section 3.1, without any further approximation. 

3.3 Higher-order schemes 

Higher-order extensions can be obtained from the above first-order algorithm by the standard tech- 
nique of composition rules. The usual symmetric composition is obtained by first taking a step of 
the above integrator with step size h/2 followed by taking its steps in reverse order. The resulting 
scheme looks as follows (hereAo = A(?o),A 1/ / 2 =A(to + |),Ai =A(tQ + h)): 

K 1/2 = U S + (A 1/2 -A )Vo, 
(U 1/2 ,S 1/2 ) = QR(K l/2 ), 
So = S\/2 - Uy 2 (A 1 / 2 -Aq)Vo, 
L^VoS^ + iAi-Ao^U,^ 

(Vi,SJ")=QR(Li), (12) 
S 1/2 = S 1 -U? /2 (A 1 -A 1/2 )V U 

K 1 = U 1/2 S 1/ 2+{A 1 -A 1/2 )Vu 
(U l ,S l )=QR(K l ), 

Fi=l/iSiVi T . 

This symmetrized splitting is a second-order scheme for (|6]). Higher-order schemes are obtained 
by suitable further compositions; see, e.g., [ffl|9]|. 

3.4 The integrator for matrix differential equations 

The basic first-order scheme extends straightforwardly to an explicit method for the low-rank ap- 
proximation of solutions A(t) of matrix differential equations A = F(A), where now A(t) is not 
known beforehand. The only change is that AA =A(ti) — A(?o) is replaced, in a way resembling 
the explicit Euler method, with 

AA = hF(Y ). 

The symmetric composition with the adjoint method now yields an implicit method. An explicit 
method of order 2 is obtained by first taking a step with the basic first-order integrator, which 
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yields an approximation Y\ , and then to take a step with the above second-order method in which 
A(t) is replaced, at t = t + Oh, with the linear function B(t + 9h) = (l- 6)F(Y ) + 0F(Y{), and 
correspondingly A(t) with the quadratic function Yq + Jf B(s) ds, that is, 

A(t + 0h)^Y + h [° B(t o + ^h)d^ = Y o + ^e(2-O)F(Y o ) + ^0 2 F(Y l ). 

JO 2 2 

Higher-order methods can again be obtained by composition of steps of the second-order method. 



4 Robustness under over-approximation 

The equations of motion Q break down when S becomes singular, and standard numerical in- 
tegrators applied to (|5]) run into problems with stability and accuracy when S is ill-conditioned. 
Such a situation arises when the matrix A(t) to be approximated by a rank-r matrix has rank less 
than r, or is close to a rank-deficient matrix. It is a remarkable property of the integrator proposed 
here that it behaves much better than a standard integrator applied to ([5]) in such a situation of 
over-approximation. On the one hand, this is already apparent from the observation that there is 
no matrix inversion in the algorithm. There is in fact more to it. 

The following result depends on the ordering of the splitting of the projector ([8]), so that we first 
compute K = US, then S, then L = VS T . For a different ordering, such as computing subsequently 
K, L, S, the following surprising exactness result is not valid. 

Theorem 4.1. Suppose that A(t) has rank at most rfor all t. With the initial value Yq = A(to), the 



splitting algorithm of Section 3.2 is then exact: Y\ =A(t\). 

Proof. We decompose A(t) = U (t)S(t)V(t) T , where both U (t) and V (t) have r orthonormal columns, 
and S(t) is an r x r matrix. We assume that V(t\) T V(to) is invertible. If this is not satisfied, then 
we make the following argument with a small perturbation of A(t\) such that V(ti) T V(to) becomes 
invertible, and let the perturbation tend to zero in the end. 

The first substep of the algorithm, starting from Yq = UoSoV Q T = U (to)S(to)V(to) T =A(to), yields 

U l S 1 =A(t 1 )V = U(t l )S(h)(V(t l ) T V(t )), 

so that the range of 

M =A{h) = (t/(? 1 )5(? 1 ))y(r 1 ) T = U l S l (V(t l ) T V(t )y l V(h) T 

is contained in the range of U\, and hence we have 

UiUiAi = A\ , as well as AqVoV t = A . 

Using the formulas of the splitting scheme we then calculate 

Fi = t/iSiVj 7 

= [/i5 y T + C/if/ 1 T AA 

= Ui Si V T - Ui uJMVqV q t + Ui U^AA 

= U S V T + AAV V T - Ui Uj AAV V T + U\ U^AA 

= Ao+AtfoVj -A -A l V V T + U 1 U?A +A 1 -U 1 U?A =A h 

which is the stated result. □ 
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Consider now a small perturbation to a matrix of rank less than r: with a small parameter £, 
assume that (with primes as notational symbols, not derivatives), 

A(t) =A'(t) + £A"(t) withrank(A / (/ 1 )) = q< r, 

where A' and A" and their derivatives are bounded independently of £. We factorize 

A' {t) = U'(t)S'(t)V'(t) T 

with U'(t) and V'(t) having q orthonormal columns and with an invertible q x q matrix S'(t). 

We apply the splitting integrator for the dynamical rank-r approximation of A(t) with starting 
value 

Y =A'(t ) + £A^ rank(y )=r, 

where Aq is bounded independently of £ (but may differ from A" (to)). We compare the result of 
the rank-r algorithm with that of the rank-g algorithm starting from 

? = A' (to) + eA'o', rank(fo) = q < r. 

Theorem 4.2. In the above situation, let Y n and Y n denote the results of n steps of the splitting 
integrator for the rank-r approximation and rank-q approximation, respectively, applied with step 
size h. Then, as long as to + nh < T, 

\\Y n -Y n \\<C(e + h), 

where C is independent ofn, h and £ (but depends onT — to). 

We note that by the standard error estimates of splitting methods, the integration error of the 
rank-g approximation is Y n — Y(t n ) = 0(h p ), uniformly in £, for the integrator of order p. Further- 
more, it follows from the over-approximation lemma in Q, Section 5.3, that the difference of the 
rank-r and rank-g approximations is bounded by Y(t) —Y(t) = 0(e). 

Proof, (a) We factorize 

Y = UoSoV T , 

where U = (U^Ug) G R mxr = R mxq x R mx ( r -9) and V = (V^Vtf) G R nxr = R nxq x R nx ( r -i) 
have orthonormal columns. Sq is chosen as an r x r matrix in block-diagonal form 

So = ^ 5 °„ j with S'q = S'(t ) and S' ' = 0(e). 

We consider the differential equation in the first splitting step, Y[(t) = A(t)P^, y r n\y We factorize 
(omitting the subscript / and the argument t) 

Y = USV T , 

where U and V have r orthonormal columns, so that we have the equation 

USV T + USV T + USV T =AVV T . (13) 
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As is shown in the proof of Lemma 5.4 of (see also B3), the decomposition becomes unique if 
we impose that S stays block-diagonal, 



S 

and 



S' 
S">' 



U T U = H, V T V = K 

with rxr matrices of the block form 

which are skew-symmetric, H12 = —#21 and K12 = —K^ ■ With the corresponding decompositions 
U = (£/', U") and V = (V', V") the proof of Lemma 5.4 of [5] yields that 

H n = -(S')- T V' T A T U" + 0{e), K n = -(S')- l U ,T AV" + 0(e), 



and in particular H\2 = 0(1), K\2 = 0(1). Multiplying (13) with U T from the left and with V from 
the right, we obtain 



HS + S + SK T = U T AV, 



which yields on the diagonal blocks 



S = U ,T AV' : s" = U" T AV". 



From ( |T3| ) we further obtain 

(u's')' + u"s"kJ 2 ='av' 

.1 

and, using also the above equation for S , 

S'V T +H n S"V" T = U' T AV"V" T . 

We show in part (b) of the proof below that U ;;T AV ;; = 0(e + h) (but we cannot conclude the same 
for U' T AV"). In summary, we get (indicating now the subscript /) 

Sj = Uj T AVi, Sj = 0(e + h) 
(Uj^)' =AV; + 0(£ + h) 
S'jv'J = Uj T AVj'Vi' T + 0(e + h). 
For the second step in the splitting we obtain similarly 

Sjj = -U^AVij, Sa = 0(e + h) 
UjjS'n = -U'ip'ffAV'n + 0(e + h) 

S'nVn = -UffAVnVF + 0(e + h) , 
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and for the third step we obtain 

S III = UlJ I AVj II , s" m = 0(e + h) 

UmS'm = U'/jjUfiAVln + 0(e + h) 

{S' ni v7 n )' = U' I ] I A + 0{e + h). 

Comparing these equations with those for S, [7, V in the rank-g splitting method, it follows that 

S\=S , 1I1 (t l )=S l +0(h£ + h 2 ) 
U[ = U' m {h) =U X + 0(he + h 2 ) 
Vl = v; il (t l )=V l +0(he + h 2 ). 

Since stable error propagation in the rank-g splitting method is obtained by the standard argument 
using the Lipschitz continuity of the right-hand side of the differential equation (uniformly in £), 
we conclude to the assertion of the theorem, 
(b) We decompose the rank-g matrix 

A'{t) = U'(t)S'(t)V'(t) T with U'(t)6'(t) = 0, V (t)9'(t) = 

and U'(to) = Uq + 0(e + h), V'(t Q ) = Vq + 0(e + h). This choice of initial factors is possible 
because of the condition A' (to) = Yq + 0(e + h), and the factorization at later t exists by solving 
the rank-g differential equations For t < t < t + h we have U'(t) = U'(t ) + 0(h) = U'(t ) + 
0(e + h) = U'(t) + 0(e + h) and V'(t) = V'(t) + 0(e + h) , and hence 

A'(t) = 6'(t)S'(t)V'(t) T + U'(t)§'(t)V'(t) T + U'(t)S'(t)9'(t) T + 0(e + h). 
Since U"(t) T U'(t) = and V"(t) T V'(t) = 0, this yields the desired bound 

U"(t) T A(t)V"(t) = 0(£ + h), 
and the proof is complete. □ 

5 Numerical experiments 
5.1 Problem setting 

The example is taken from 0. We generate time-dependent matrices A (t) as 

A(t) = Qi(t)(A l +exp(t)A 2 )Q 2 (t), 

where the orthogonal matrices Qi(t) and Q 2 (t) are generated as the solutions of differential equa- 
tions 

Qi = TiQu i=l,2 

with random skew-symmetric matrices 7]. The matrices Ai,A 2 are generated as follows. First, 
we generate a 10 x 10 matrix as an identity matrix plus a matrix with random entries distributed 
uniformly over [0,0.5]. This matrix is then set as a leading block of a 100 x 100 matrix. After 
that, we add a perturbation to this enlarged matrix. The perturbation is generated as a matrix with 
uniformly distributed entries on [0,e]. For small £ this matrix is close to a rank- 10 matrix. If the 
rank of the approximation is chosen larger than 10, the norm of 5 1 in ^ will be large and this 
leads to instability with standard time discretizations. 
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5.2 Numerical comparisons 

We will compare the following schemes. The first scheme is the standard implicit midpoint rule 
combined with fixed point iteration applied directly to the system ([5]). We test the proposed splitting 
schemes with different orders of splitting and with/without symmetrization. We denote by KSL the 
scheme of Section 3.2, where first K = US is updated, then S, then L = VS T . By KLS we denote 
the scheme where first K, then L, then S are updated. We thus consider the following methods: 

1 . KLS scheme (first order) 

2. KLS scheme with symmetrization (second order) 

3. KSL scheme (first order) 

4. KSL scheme with symmetrization (second order) 

We are interested in the approximation errors \\Y(t) — A(t)\\ for each particular scheme and differ- 
ent values of r and £. The results are presented in Figure [T] where we also plot the error of the best 
rank-r approximation to A{t) computed by SVD. Note that both KSL schemes perform remarkably 
better in the overapproximation case (subplot d). The midpoint rule is unstable in case d), whereas 
both KLS schemes have significantly higher error than the KSL schemes. All computations are 
done with constant step size h = 10 -3 . 

To test the convergence properties of the schemes with respect to the timestep h, we have com- 
puted the numerical order of the different schemes using the Runge rule: 

\\y(h) -y(h/2) || / \\y(h/2) - y(h/4)\\ w 2 P in case of order p. 





P 


Appr. err. 




P 


Appr. err. 


Midpoint 


2.0023 


0.2200 


Midpoint 


2.0024 


0.0188 


KLS 


1.0307 


1.8133 


KLS 


1.0309 


1.8030 


KLS(symm) 


1.8226 


0.2215 


KLS(symm) 


1.8231 


0.0324 


KSL 


1.0089 


0.2188 


KSL 


1.0082 


0.0002 


KSL(symm) 


2.005 


0.2195 


KSL(symm) 


2.0049 


0.0002 


Table 1: £ = 10 3 , 


r= 10 


Table 


2: e = 10 


- 6 ,r= 10 




P 


Appr. err. 




P 


Appr. err. 


Midpoint 


0.0001 


0.1006 


Midpoint 




failed 


KLS 


0.8154 


1.4224 


KLS 


0.9633 


1.3435 


KLS(symm) 


1.4911 


0.3142 


KLS(symm) 


0.3127 


1.5479 


KSL 


1.0354 


0.0913 


KSL 


1.0362 


9.1316e-05 


KSL(symm) 


1.9929 


0.0913 


KSL(symm) 


1.993 


9.1283e-05 


Table 3: e = 10 3 , 


r = 20 


Table 


4: e = 10 


~ 6 ,r = 20 



The approximation error listed is the error at t = 1 with respect to the given matrix A(t), not with 
respect to the solution Y(t) of rank r of the differential equations ([5]). In the overapproximation 
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*- -* Best approx. 

• — • Midpt. rule 

- - KLS 

♦ — ♦ KLS(symm) 

>^ KSL 

KSL(symm) 





0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 

t t 



c)e = l(T 3 ,r = 20 d)£ = 10 _6 ,r = 20 



Figure 1 : Dynamical low-rank approximation error for different schemes and different values of 
approximation rank r and perturbation size £. 
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case both KLS schemes perform much worse, and the symmetrized version loses its second order. 
The KSL scheme, that is, the scheme of Section 3.2, and its symmetrized version (see Section 3.3) 
clearly outperform the other methods. 

The last test describes the stability of the schemes with respect to the time step h. In Figure [2] 
we plot the approximation error at time t = 1 for both KSL schemes and the midpoint rule for 
different h ranging from 10 1 to 10~ 3 . The rank r was set to 20 (overapproximation case) and the 
noise level was chosen to be e = 10~ 3 . The midpoint rule becomes unstable for large values of h, 

xlO" 2 



9.05 

o 

CD 

.2 9.00 



I 8.95 
< 

8.90 

10" 3 10" 2 10 _1 

h 




Figure 2: Dynamical low-rank approximation error for different schemes versus stepsize h, for 
r = 20 and £ = 10~ 3 

whereas both KSL schemes give a good result for the whole range of stepsizes. It is interesting to 
compare this value with the best low -rank approximation of A(l). The error of the best rank- 10 
approximation of A(l) (computed by the SVD) is approximately 10 and the error of the best 
rank-20 approximation is approximately 3 ■ 10~ 3 . The dynamical low-rank approximation thus 
captures the "smooth" component of the solution. 

6 Conclusion and perspectives 

We have proposed and analyzed a fully explicit, computationally inexpensive integrator for the dy- 
namical low-rank approximation of time-dependent matrices that is based on splitting the projector 
onto the tangent space of the low-rank manifold. The integrator has remarkable robustness under 
over-approximation with a too high rank. While standard explicit and implicit integrators break 
down in such a situation, the integrator proposed here does not suffer from the ill-conditioning or 
singularity of the small matrix factor in the orthogonal low-rank factorization. 

The robustness under overapproximation enables one to control the rank adaptively. Lowering 
the rank is trivial, but thanks to the robustness under a reduced rank we are able to raise the rank 
in a natural way, continuing the computations with the higher rank starting from the values of 
lower rank. This is not possible with standard integrators applied to the higher-rank differential 
equations, which would have to start from a singular matrix factor, in which case the differential 
equations are not well-defined. 

Another application of the proposed integrator is in optimization algorithms on a low-rank man- 
ifold, such as eg or Newton's method. There, an update A + AA to a low-rank iterate A has to 
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be truncated (or retracted in another terminology) back to the given low rank. This can be done 
efficiently by applying one step of our integrator to A + tAA at t = 1 . 

The integrator can be extended to the low-rank approximation of tensors in the tensor train and 
hierarchical Tucker formats. This extension will be studied in forthcoming work. 
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